/***************************************************************************
 * matrices.c  -  description
 * ----------------------------
 * begin                : unknown
 * author               : Boyan Bejanov
 * email                : bbejanov@ualberta.ca
 ***************************************************************************/

/****************************************************************************
 *  Matrices functionality routines definitions                             *
 ****************************************************************************/

#include "matrices.h"
#include "errlog.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>


#if defined(USE_HYPRE)
#include <HYPRE.h>
#endif

#if defined(USE_SPSOLVER)
#include <sunperf.h>
#endif

/****************************************************************************
 * MATRIX and VECTOR memory management routines
 ****************************************************************************/

VECTOR allocV(int dim)
/****************************************************************************
 * Routine:     allocV
 * ------------------------------------------------------------------------
 * Purpose:
 *      Allocates memory for a VECTOR and initializes it to 0
 * Parameters:
 *      int dim     - dimension of VECTOR
 * Return value:
 *      Pointer to the newly allocated vector
 * Algorithm:
 *      calloc is called to allocate memory and set it to 0.
 *      if calloc fails program is terminated with error
 ****************************************************************************/
{
    VECTOR ret=(VECTOR)calloc(dim, sizeof(double));
    if(ret==NULL)
        error_msg("cannot allocate vector.");
    return ret;
}

void freeV(VECTOR *V)
/****************************************************************************
 * Routine:     freeV
 * ------------------------------------------------------------------------
 * Purpose:
 *      Releases memory used by a vector and sets VECTOR to NULL
 * Parameters:
 *      VECTOR *V   - pointer to the vector to be destroyed
 * Algorithm:
 *      call free to release memory and then set V to 0(NULL)
 ****************************************************************************/
{
    free(*V);
    (*V)=NULL;
}

MATRIX allocM(int rows, int cols)
/****************************************************************************
 * Routine:     allocM
 * ------------------------------------------------------------------------
 * Purpose:
 *      Allocates memory for a MATRIX and initializes it to 0
 * Parameters:
 *      int rows, cols  - dimensions of MATRIX
 * Return value:
 *      Newly allocated matrix
 * Local variables:
 *      MATRIX M    - the matrix
 *      int i       - loop counter
 * Algorithm:
 *      First allocate array of VECTORs of size rows. These will be rows.
 *      Call allocV to allocate all rows*cols entries of the matrix.
 ****************************************************************************/
{
    MATRIX M;
    int i;
    M=(MATRIX)malloc(rows*sizeof(VECTOR));
    if(M==NULL)
        error_msg("Cannot allocate matrix.");
    M[0] = allocV(rows * cols);
    for(i=1;i<rows;i++)
        M[i]=M[i-1] + cols;
    return M;
}

void freeM(MATRIX *M)
/****************************************************************************
 * Routine:     freeM
 * ------------------------------------------------------------------------
 * Purpose:
 *      Releases memory used by MATRIX.
 * Parameters:
 *      MATRIX *M   - pointer to MATRIX to be destroyed
 * Algorithm:
 *      Release the memory used by the elements. Release the vector of rows
 ****************************************************************************/
{
    free((*M)[0]);
    free(*M);
    (*M)=NULL;
}

void CopyVector(const double *S, VECTOR T, int Length)
/****************************************************************************
 *  Routine:     CopyVector
 * ------------------------------------------------------------------------
 *  Purpose:
 *      Copies the contents of one VECTOR into another
 *  Parameters:
 *      VECTOR S    - source
 *      VECTOR T    - target
 *      int Length  - dimension of vector (# of elements to be copied)
 ****************************************************************************/
{
    memcpy(T, S, Length*sizeof(double));
    return;
}

void CleanVector(VECTOR V, int Length)
/****************************************************************************
 *  Routine:     CleanVector
 * ------------------------------------------------------------------------
 *  Purpose:
 *      Sets the contents of a VECTOR to 0
 *  Algorithm:
 *      simply call memset to do the job
 ****************************************************************************/
{
    memset(V, 0, Length*sizeof(double));
    return;
}


/****************************************************************************
 * Sparse_Row
 ****************************************************************************/

/*  Structures of the following type will contain information describing one */
/*  row of a sparse matrix. This includes an array of the column indices of  */
/*  the nonzero elements in the row, an array of the corresponding values  */
/*  and the number of entries already allocated */
typedef struct _SPARSE_ROW
{
    int EntriesCount;    /* Number of nonzero matrix entries in this row */
    int *Columns;        /* Array of column indices of nonzero elements */
    double *Values;          /* Array of the values of these entries */
} Sparse_Row;

static Sparse_Row *AllocSparse_Rows(int rows, int sym)
{
    int i;
    Sparse_Row *Row, *R;
    /* All rows are by default empty, i.e. their EntriesCount's are 0
     * ( as well as the pointers Columns and Values ).
     * Allocate memory with calloc, so it is initialized with 0 */
    R=Row=(Sparse_Row*)calloc(rows, sizeof(Sparse_Row));
    if(Row==NULL)
        error_msg("Cannot allocate sparse matrix.");
    /* In Symmetric matrices we put the diagonal element at the top of the
     * array, because it is treated in a different way. This way we don't
     * have to check if the current element is diagonal -- we just handle
     * the first entry as diagonal and the rest of the elements in the row
     * as off-diagonal. Here we allocate the diagonal elements: */
    if(sym)
    {
        /* Row=(Sparse_Row*)M->SPARSE_DATA; */
        for(i=0;i<rows;i++, Row++)
        {
            if( ((Row->Columns=(int*)malloc(sizeof(int)))==NULL)
            || ((Row->Values=(double*)malloc(sizeof(double)))==NULL)
            )
                error_msg("cannot allocate memory for sparse matrix.");
            Row->EntriesCount=1;
            Row->Columns[0]=i;
            Row->Values[0]=0.0;
        }
    }
    return R;
}

static void CleanSparse_Row(Sparse_Row *R, int rows)
{
    int i;
    for(i=0;i<rows;i++, R++)
        memset(R->Values, 0, R->EntriesCount * sizeof(double));
}

static void DestroySparse_Row(Sparse_Row *Row, int rows)
{
    int i;
    Sparse_Row *R = Row;
    if(Row!=NULL)
    {
        for(i=0;i<rows;i++, Row++)
        {
            if(Row->EntriesCount)
            { free(Row->Columns); free(Row->Values); }
        }
        free(R);
    }
    
}

static void PrintSRow(FILE *F, Sparse_Row *R, int rows)
{
    int i, j, *c;
    double *v;
    for(i=0;i<rows;i++, R++)
    {
        c = R->Columns;
        v = R->Values;
        for(j=0;j<R->EntriesCount;j++, v++, c++)
            fprintf(F, "%d,%d,%.20g\n", i, (int)*c, (double)*v);
    }
}

static void _SM_IncreaseRowArrays(Sparse_Row *R, int incr)
/****************************************************************************
 * Routine:     _SM_IncreaseRowArrays
 * ------------------------------------------------------------------------
 * Purpose:
 *      Increases the amount of memory allocated to the arrays in a row
 *      NOTE: This routine is internal for this module /not visible to the
 *      user/ and is only called by SM_FindEntry with incr=1 (see below).
 * Parameters:
 *      Sparse_Row *R   - pointer to the row whose arrays are to be altered
 *      int inct        - increment in number of matrix entries
 * Local variables:
 *      void *ptr       - pointer. Address of newly allocated memory block
 *              will be assigned to this variable, checked, and then
 *              assigned to its corresponding field in row structure.
 *      int ns          - the new size of the arrays
 * Algorithm:
 *      Call realloc twice (for Columns and Values) to do the job.
 ****************************************************************************/
{
    void *ptr;
    int ns;
    
    ns = R->EntriesCount + incr;
    ptr = realloc(R->Columns, ns*sizeof(int));
    if( ptr == NULL )
        error_msg("cannot allocate memory for sparse matrix.");
    R->Columns=(int*)ptr;
    ptr = realloc(R->Values, ns*sizeof(double));
    if( ptr == NULL )
        error_msg("cannot allocate memory for sparse matrix.");
    R->Values=(double*)ptr;
    R->EntriesCount = ns;
}

static double* _SM_FindEntrySparse_Row(Sparse_Row *R,
        int row, int col, int sym, int add)
{
    int i;
    if( sym && (row > col) )
    {
        i = row;
        row = col;
        col = i;
    }
    
    R+=row;
    for(i=0;i<R->EntriesCount;i++)
        if(R->Columns[i]==col)
            return &R->Values[i];
    if(!add)
        return NULL;
    _SM_IncreaseRowArrays(R, 1);
    /*    i=R->EntriesCount-1; */
    R->Columns[i]=col;
    R->Values[i]=0.0;
    return &R->Values[i];
    
}


static void _SM_MultSymmV_Sparse_Row(Sparse_Row *Rows,
        int dim, const double *vect, double **res)
        /****************************************************************************
         * Routine:     _SM_MultSymmV
         * ------------------------------------------------------------------------
         * Purpose:
         *      Multiplies a symmetric Sparse_Matrix by a vector
         * Parameters:
         *      Sparse_Row *Rows    - the rows data of the sparse matrix
         *      int dim          - the dimesion of the result vector
         *      VECTOR vect         - the vector to be multiplied
         *      VECTOR *res         - the result
         * Return value:
         * Local variables:
         *      int i, j         - loop counters
         *          *Cols           - the array of column indices in a row
         *      double *Vals         -  the array of values in a row
         * Algorithm:
         *      ...
         ****************************************************************************/
{
    int i, j, *Cols;
    double *Vals, *r;
    const double *v;
    Sparse_Row *R=Rows;
    memset(*res, 0, dim*sizeof(double));
    r = *res, v = vect;
    for(i=0;i<dim;i++, R++, r++, v++)
    {
        Vals=R->Values;
        Cols=R->Columns;
        (*r)+=(*v)*(*Vals);
        for(j=1, Cols++, Vals++;j<R->EntriesCount;j++, Cols++, Vals++)
        {
            (*r)+=vect[*Cols]*(*Vals);
            (*res)[*Cols]+=(*v)*(*Vals);
        }
    }
}

static void _SM_MultV_Sparse_Row(Sparse_Row *Rows,
        int dim, const double *vect, double **res)
        /****************************************************************************
         * Routine:     _SM_MultV
         * ------------------------------------------------------------------------
         * Purpose:
         *      Multiplies a non-symmetric Sparse_Matrix by a vector
         * Parameters:
         *      Sparse_Row *Rows    - the rows data of the sparse matrix
         *      int dim          - the dimesion of the result vector
         *      VECTOR vect         - the vector to be multiplied
         *      VECTOR *res         - the result
         * Return value:
         * Local variables:
         *      int i, j         - loop counters
         *          *Cols           - the array of column indices in a row
         *      double *Vals         -  the array of values in a row
         * Algorithm:
         *      ...
         ****************************************************************************/
{
    int i, j, *Cols;
    double *Vals;
    Sparse_Row *R=Rows;
    memset(*res, 0, dim*sizeof(double));
    for(i=0;i<dim;i++, R++)
    {
        Cols=R->Columns;
        Vals=R->Values;
        for(j=0;j<R->EntriesCount;j++, Cols++)
        {
            (*res)[i]+=vect[*Cols]*Vals[j];
        }
    }
}

static void _SM_MultTransV_Sparse_Row(Sparse_Row *Rows,
        int rows, int cols, const double *vect, double **res)
        /****************************************************************************
         * Routine:     _SM_MultTransV
         * ------------------------------------------------------------------------
         * Purpose:
         *      Multiplies the transpose of a non-symmetric matrix by a vector
         * Parameters:
         *      Sparse_Row *Rows    - the rows data of the sparse matrix
         *      int dim          - the dimesion of the result vector
         *      VECTOR vect         - the vector to be multiplied
         *      VECTOR *res         - the result
         * Return value:
         * Local variables:
         *      int i, j         - loop counters
         *          *Cols           - the array of column indices in a row
         *      double *Vals         -  the array of values in a row
         * Algorithm:
         *      ...
         ****************************************************************************/
{
    int i, j, *Cols;
    double *Vals;
    Sparse_Row *R=Rows;
    memset(*res, 0, cols*sizeof(double));
    for(i=0;i<rows;i++, R++)
    {
        Cols=R->Columns;
        Vals=R->Values;
        for(j=0;j<R->EntriesCount;j++, Cols++)
        {
            (*res)[*Cols]+=vect[i]*Vals[j];
        }
    }
}

/*****************************************************************************/

/****************************************************************************
 * Sparse_CSC  - compressed sparse column format
 ****************************************************************************/

/* NOTE: http://docs.sun.com/source/816-2463/plug_matrices.html#pgfId-12370 */

/****** IMPORTANT NOTE ************ : all indices in Sparse_CSC are 1-based,
 *  since they must be compatible with fortran routines
 ****************************************************************************/


#define f2c(ind)        ((long)((long)(ind)-1))
#define c2f(ind)        ((long)((long)(ind)+1))

#define SPSOLVER_PRINT_CONDITION_NUMBER
/*
 #define CSC_SORT_COLS
 #define f2c(ind)        (ind)
 #define c2f(ind)        (ind)
 */

/* possible values are "nat", "mmd", "gnd", and "usr" */
static char ordmth[] = "mmd";  

typedef struct _SPARSE_CSC
{
    char    mtrtyp[3];   /* one of "sp", "ss", "su", or "uu" */
    long    neqns;
    long    *colstr;
    long    *rowind;
    double  *values;
    double  *handle;
    long    **aux_data;
    struct
    {
        unsigned int   structure : 1;
        unsigned int   assemble  : 1;
        unsigned int   changed   : 1;
    } flags;
} Sparse_CSC;

static void BeginAssembleSparse_CSC(Sparse_CSC *data)
{
    data->flags.structure = 0;
    data->flags.assemble = 1;
    data->flags.changed = 0;
}

static void EndAssembleSparse_CSC(Sparse_CSC *data)
{
    if( data->flags.changed == 1)
    {
        warning_msg("%s: structure of matrix changed during reassembling.\n", __func__);
        #if defined(USE_SPSOLVER)
        if(data->handle)
        {
            long ierr = 0l;
            dgssor_64_(ordmth, data->handle, &ierr, 3);
            if(ierr != 0) error_msg("%s: ierr = %d", "dgssor_64_", ierr);
        }
        #endif
    }
    #if defined(USE_SPSOLVER)
    if(data->handle)
    {
        long ierr = 0;
        dgssfa_64_(&data->neqns, data->colstr, data->rowind, data->values,
                data->handle, &ierr);
        if(ierr != 0) error_msg("%s: ierr = %d", "dgssfa_64_", ierr);
    }
    #endif
            ;
    data->flags.changed = 0;
    data->flags.structure = 0;
    data->flags.assemble = 0;
}

static void BeginStructureSparse_CSC(Sparse_CSC *data)
{
    data->flags.structure = 1;
    data->flags.assemble = 0;
    data->flags.changed = 0;
    
    #if defined(USE_SPSOLVER)
    if(data->handle)
    {
        long ierr = 0L;
        dgssda_64_(data->handle, &ierr);
        free(data->handle);
        data->handle = NULL;
    }
    #endif
            ;
    if((data->aux_data=calloc(data->neqns, sizeof(long *)))==NULL)
        goto alloc_error;
    
    for(long i=0; i<data->neqns; i++)
    {
        data->colstr[i] = 1;
        if((data->aux_data[i]=calloc(1, sizeof(long)))==NULL)
            goto alloc_error;
        data->aux_data[i][0]=c2f(i);
    }
    data->colstr[data->neqns] = 0;
    
    return;
    alloc_error:
        error_msg("%s: memory allocation error.\n", __func__);
}

static void EndStructureSparse_CSC(Sparse_CSC *data)
{
    long i;
    data->flags.structure = 0;
    data->flags.assemble = 0;
    data->flags.changed = 0;
    long total = data->colstr[0];
    data->colstr[0] = c2f(0);
    for(i=1; i<data->neqns; i++)
    {
        long save = data->colstr[i];
        data->colstr[i] = c2f(total);
        total += save;
    }
    data->colstr[data->neqns] = c2f(total);
    
    void *ptr;
    ptr=realloc(data->values, total*sizeof(double));
    if(ptr==NULL) goto alloc_error;
    data->values=ptr;
    
    ptr=realloc(data->rowind, total*sizeof(long));
    if(ptr==NULL) goto alloc_error;
    data->rowind=ptr;
    
    for(i=0;i<data->neqns;i++)
    {
        memcpy(data->rowind+f2c(data->colstr[i]), data->aux_data[i],
                (data->colstr[i+1]-data->colstr[i])*sizeof(long));
        free(data->aux_data[i]);
    }
    free(data->aux_data);
    data->aux_data=NULL;
    
    data->flags.structure = 0;
    data->flags.assemble = 0;
    data->flags.changed = 0;
    
    #if defined(USE_SPSOLVER)
    assert(data->handle == NULL);
    data->handle = calloc(150, sizeof(double));
    if(data->handle==NULL) goto alloc_error;
    
    long ierr = 0L, outunt = 6, msglvl = 0;
    dgssin_64_(data->mtrtyp, "n", &data->neqns, data->colstr, data->rowind,
            &outunt, &msglvl, data->handle, &ierr, 2, 1);
    if(ierr != 0) error_msg("%s: ierr = %d", "dgssin_64_", ierr);
    
    dgssor_64_(ordmth, data->handle, &ierr, 3);
    if(ierr != 0) error_msg("%s: ierr = %d", "dgssor_64_", ierr);
    #endif
            
            return;
    alloc_error:
        error_msg("%s: cannot allocate memory for sparse matrix.", __func__);
}




static Sparse_CSC *AllocSparse_CSC(long neqns, int sym, char *mtrtyp)
{
    /* here sym refers to matrix values.  If matrix is structurally symmetric,*/
    /* but unsymmetric values, then sym must be FALSE */
    Sparse_CSC *data;
    data = calloc(1, sizeof(Sparse_CSC));
    if(data==NULL)
        error_msg("Cannot allocate Sparse_Matrix in CSC storage.\n");
    data->neqns = neqns;
    data->mtrtyp[0]=tolower(mtrtyp[0]);
    data->mtrtyp[1]=tolower(mtrtyp[1]);
    data->mtrtyp[2]=0;
    if( strcmp(data->mtrtyp, "sp")!=0 &&
            strcmp(data->mtrtyp, "ss")!=0 &&
            strcmp(data->mtrtyp, "su")!=0 &&
            strcmp(data->mtrtyp, "uu")!=0 )
        error_msg("Invalid CSC matrix type %s\n", data->mtrtyp);
    if( sym && (strcmp(data->mtrtyp, "sp")!=0 && strcmp(data->mtrtyp, "ss")!=0) )
        error_msg("%s: symmetric matrix must have symmetric storage.\n", __func__);
    if( (!sym) && (strcmp(data->mtrtyp, "sp")==0 || strcmp(data->mtrtyp, "ss")==0) )
        error_msg("%s: unsymmetric matrix cannot have symmetric storage.\n", __func__);
    data->colstr = malloc((1+neqns)*sizeof(long));
    /* allocate diagonal entries */
    data->rowind = malloc(neqns*sizeof(long));
    data->values = calloc(neqns, sizeof(double));
    data->handle = NULL;
    if(data->colstr==NULL || data->rowind==NULL || data->values==NULL)
        error_msg("Cannot allocate Sparse_Matrix in CSC storage.\n");
    for(long i=0;i<neqns;i++)
    {
        data->colstr[i]=c2f(i);
        data->rowind[i]=c2f(i);
    }
    data->colstr[neqns]=c2f(neqns);
    return data;
}

static void CleanSparse_CSC(Sparse_CSC *data)
{
    memset(data->values, 0, sizeof(double) * f2c(data->colstr[data->neqns]) );
}

static void DestroySparse_CSC(Sparse_CSC *data)
{
    if(data)
    {
        if(data->colstr) free(data->colstr);
        if(data->rowind) free(data->rowind);
        if(data->values) free(data->values);
        if(data->handle)
        {
            #if defined(USE_SPSOLVER)
            long ierr = 0L;
            dgssda_64_(data->handle, &ierr);
            #endif
                    ;
            free(data->handle);
        }
        free(data);
    }
}

static void PrintSparse_CSC(FILE *F, Sparse_CSC *data)
{
    long col, i;
    for(col=0;col<data->neqns;col++)
        for(i=f2c(data->colstr[col]);i<f2c(data->colstr[col+1]);i++)
            fprintf(F, "%ld,%ld,%.20g\n",
                    f2c(data->rowind[i]), col, data->values[i]);
}

/*
 * static void _SM_IncreaseRowArrays(Sparse_Row *R, int incr)
 * {
 * void *ptr;
 * int ns;
 *
 * ns = R->EntriesCount + incr;
 * ptr = realloc(R->Columns, ns*sizeof(int));
 * if( ptr == NULL )
 * error_msg("cannot allocate memory for sparse matrix.");
 * R->Columns=(int*)ptr;
 * ptr = realloc(R->Values, ns*sizeof(double));
 * if( ptr == NULL )
 * error_msg("cannot allocate memory for sparse matrix.");
 * R->Values=(double*)ptr;
 * R->EntriesCount = ns;
 * }
 */

static double* _SM_FindEntrySparse_CSC(Sparse_CSC *data,
        long row, long col, int add)
{
    long i;
    void *ptr;
    
    if(data->mtrtyp[0]=='s' && (data->mtrtyp[1]=='p' || data->mtrtyp[1]=='s'))
    {
        /* storage is symmetric and values are symmetric */
        /* we only store the lower triangular part, so ...  */
        if(row < col)
        { /* swap row <-> col */ i = row; row = col; col = i; }
    }
    
    if(data->flags.structure)
    {
        for(i=0; i<data->colstr[col]; i++)
            if(row == f2c(data->aux_data[col][i]))
                return data->values;
        long foo = ++data->colstr[col];
        ptr = realloc(data->aux_data[col], foo*sizeof(long));
        if(ptr==NULL) goto alloc_error;
        data->aux_data[col] = ptr;
        #ifdef  CSC_SORT_COLS
                for(i=foo-1;i>0;i--)
                {
            if(row < f2c(data->aux_data[col][i-1]))
                data->aux_data[col][i] = data->aux_data[col][i-1];
            else
            {
                data->aux_data[col][i] = c2f(row);
                break;
            }
                }
        if(i==0) data->aux_data[col][0] = c2f(row);
        #else
                data->aux_data[col][foo-1] = c2f(row);
        #endif
                return data->values;
    }
    
    
    for(i=f2c(data->colstr[col]); i<f2c(data->colstr[col+1]); i++)
    {
        if(row == f2c(data->rowind[i]))
            return data->values + i;
        #ifdef CSC_SORT_COLS
                if(row < f2c(data->rowind[i]))
                    break;
        #endif
    }
    
    if(!add) return NULL;
    
    data->flags.changed = 1;
    
    /* add entry (row,col) at position i in arrays, push next entries back */
    for(col++; col<=data->neqns; col++) data->colstr[col]++;
    
    long nentries = f2c(data->colstr[data->neqns]);
    
    ptr=realloc(data->values, nentries*sizeof(double));
    if(ptr==NULL) goto alloc_error;
    data->values=ptr;
    
    ptr=realloc(data->rowind, nentries*sizeof(long));
    if(ptr==NULL) goto alloc_error;
    data->rowind=ptr;
    
    memmove(data->values+i+1, data->values+i, (nentries-i-1)*sizeof(double));
    memmove(data->rowind+i+1, data->rowind+i, (nentries-i-1)*sizeof(long));
    
    data->rowind[i]=c2f(row);
    data->values[i]=0.0;
    return data->values+i;
    
    alloc_error:
        error_msg("%s: cannot allocate memory for sparse matrix.", __func__);
}

static void _SM_MultV_Sparse_CSC(Sparse_CSC *data,
        const double *vect, double **res)
{
    long nentries = f2c(data->colstr[data->neqns]);
    long i, row, col;
    double val;
    col=0;
    memset(*res, 0, data->neqns*sizeof(double));
    if(((strcmp(data->mtrtyp, "ss")==0) || (strcmp(data->mtrtyp, "sp")==0)))
        for(i=0;i<nentries;i++)
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[row] += val * vect[col];
        if(row!=col) (*res)[col] += val * vect[row];
        }
    else
        for(i=0;i<nentries;i++)
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[row] += val * vect[col];
        }
}

static void _SM_MultTransV_Sparse_CSC(Sparse_CSC *data,
        const double *vect, double **res)
{
    long nentries = f2c(data->colstr[data->neqns]);
    long i, row, col;
    double val;
    col=0;
    memset(*res, 0, data->neqns*sizeof(double));
    if(((strcmp(data->mtrtyp, "ss")==0) || (strcmp(data->mtrtyp, "sp")==0)))
        for(i=0;i<nentries;i++)     /* symmetric is the same as straight */
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[row] += val * vect[col];
        if(row!=col) (*res)[col] += val * vect[row];
        }
    else
        for(i=0;i<nentries;i++)     /* transpose multiplication */
        {
        while(i>=f2c(data->colstr[col+1])) col++;
        row=f2c(data->rowind[i]);
        val=data->values[i];
        (*res)[col] += val * vect[row];
        }
}

/*****************************************************************************/



/****************************************************************************
 * Sparse Matrices
 ****************************************************************************/
void SM_BeginStructure(Sparse_Matrix *M)
{
    if(M->data_type[0]==0)
        return;
    else
        BeginStructureSparse_CSC(M->SPARSE_DATA);
    return;
}

void SM_EndStructure(Sparse_Matrix *M)
{
    if(M->data_type[0]==0)
        return;
    else
        EndStructureSparse_CSC(M->SPARSE_DATA);
    return;
}

void SM_BeginAssemble(Sparse_Matrix *M)
{
    CleanSparseMatrix(M);
    if(M->data_type[0]==0)
        return;
    else
        BeginAssembleSparse_CSC(M->SPARSE_DATA);
    return;
}

void SM_EndAssemble(Sparse_Matrix *M)
{
    if(M->data_type[0]==0)
        return;
    else
        EndAssembleSparse_CSC(M->SPARSE_DATA);
    return;
}


void CleanSparseMatrix(Sparse_Matrix *M)
/****************************************************************************
 *  Routine:     CleanSparseMatrix
 * ------------------------------------------------------------------------
 *  Purpose:
 *      Sets the entries of a Sparse_Matrix to 0
 *  Algorithm:
 *      Does not release any memory allocated to rows
 ****************************************************************************/
{
    if((M==NULL)||(M->rows==0)||(M->cols==0)||(M->SPARSE_DATA==NULL))
        return;
    if(M->data_type[0]==0)
        CleanSparse_Row(M->SPARSE_DATA, M->rows);
    else
        CleanSparse_CSC(M->SPARSE_DATA);
    return;
}

int CreateSM(Sparse_Matrix *M, int rows, int cols, int sym,
        char *type)
        /****************************************************************************
         * Routine:     CreateSM
         * ------------------------------------------------------------------------
         * Purpose:
         *      Allocates and initializes the necesary data structures for a
         *      sparse matrix
         * Parameters:
         *      Sparse_Matrix *M   - Receives a pointer to the matrix to be created
         *      int rows, cols     - dimensions of the new matrix
         *      int sym            - TRUE(!=0) if matrix is to be symmetric
         *                           FALSE(0) otherwise
         * Return value:
         *      TRUE if successful, or FALSE if error occured
         * Local variables:
         *      Sparse_Row *Row    -
         *      int i              - loop counters
         * Algorithm:
         *      First, assign the values of the paramenters to their corresponding
         *      fields in the sparse matrix structure. Note that a matrix is only
         *      set to be symmetric if it is square (rows==cols). If this is not
         *      true the matrix will not be symmetric even if user requests this.
         *      Next, allocate memory for the array of rows and set it to 0 by
         *      calling 'calloc'. This is important -- see below.
         *      Finally, if matrix is symmetric, allocate an entry for the diagonal
         *      element in each row
         ****************************************************************************/
{
    void *data;
    
    if ( M == NULL )
        return FALSE;
    M->rows=rows;
    M->cols=cols;
    if((rows==cols)&&sym)
        M->Symmetric=TRUE;
    else
        M->Symmetric=FALSE;
    if(type==NULL || (rows!=cols))
    {
        M->data_type[0] = 0;
        data=AllocSparse_Rows(rows, M->Symmetric);
    }
    else
    {
        M->data_type[0]=type[0];
        M->data_type[1]=type[1];
        data=AllocSparse_CSC(rows, M->Symmetric, type);
    }
    M->SPARSE_DATA = data;
    return TRUE;
}

void DestroySM(Sparse_Matrix *M)
/****************************************************************************
 * Routine:     DestroySM
 * ------------------------------------------------------------------------
 * Purpose:
 *      Releases memory used by sparse matrix
 * Parameters:
 *      Sparse_Matrix *M    - matrix to be destayed
 * Local variables:
 *      Sparse_Row *Row         -
 *      int i                   - loop counters
 * Algorithm:
 *      In a loop for each row free memory used by Columns and Values
 *      arrays. Then free memory used by the rows, and clear the matrix.
 ****************************************************************************/
{
    if ( M == NULL )
        return;
    if(M->data_type[0]==0)
        DestroySparse_Row(M->SPARSE_DATA, M->rows);
    else
        DestroySparse_CSC(M->SPARSE_DATA);
    memset(M, 0, sizeof(Sparse_Matrix));
}


static double* SM_FindEntry(Sparse_Matrix *M, int row, int col,
        int add)
        /****************************************************************************
         * Routine:     _SM_FindEntry
         * ------------------------------------------------------------------------
         * Purpose:
         *      Locates the memory address of an entry of a sparse matrix
         * Parameters:
         *      Sparse_Matrix M     - Sparse matrix
         *      int row, col        - coordinates of entry
         *      int add             - indicates wether a new entry shall be created
         *                  if the one does not exist (meaning its value is 0)
         * Return value:
         *      double*  pointer to the entry or NULL if entry is not found and
         *      add is 0
         * Local variables:
         *      Sparse_Row *R       - pointer to the row of the element
         *      int i               - loop conter
         * Algorithm:
         *      Firstly, check if dimensions are correct. Return NULL if not --
         *      this will indicate that the entry was not found in memory, which
         *      can be because entry is 0 and therefore not allocated, or because
         *      (as in this case) entry does not exist (error). If matrix is
         *      symmetric then the row must be the smaller index. Get the row of
         *      the element in R and scan its Columns array. If the entry doesn't
         *      exist then we create a new one (in case of add) or return NULL (!add)
         ****************************************************************************/
{
    if( (row>=M->rows) || (col>=M->cols) )
        error_msg("FindEntry: indices out of range\n");
    if(M->data_type[0]==0)
        return _SM_FindEntrySparse_Row(M->SPARSE_DATA, row, col,
                M->Symmetric, add);
    else
        return _SM_FindEntrySparse_CSC(M->SPARSE_DATA, row, col, add);
}

void SM_SetEntry(Sparse_Matrix *M, int row, int col, double AValue)
/****************************************************************************
 * Routine:     SM_SetEntry
 * ------------------------------------------------------------------------
 * Purpose:
 *      Assignes new value to an entry in sparse matrix
 * Parameters:
 *      Sparse_Matrix M     - the sparse matrix
 *      int row, col        - coordinates of element in matrix
 *      double AValue        - new value
 * Local variables:
 *      double *E        - pointer to the location in memory where the value
 *              of this entry is stored
 * Algorithm:
 *      Call _SM_FindEntry to locate the memory address of the entry. It is
 *      called with _add_=TRUE (see above) so new entry will be allocated if
 *      necessary.
 *      Perform the operation using the pointer returned by _SM_FindEntry.
 ****************************************************************************/
{
    double *E;
    E=SM_FindEntry(M, row, col, TRUE);
    (*E)=AValue;
    return;
}

void SM_AddToEntry(Sparse_Matrix *M, int row, int col, double AValue)
/****************************************************************************
 * Routine:     SM_AddToEntry
 * ------------------------------------------------------------------------
 * Purpose:
 *      Increases the value of an entry in sparse matrix with given value
 *      Useful for assembling
 * Parameters:
 *      Sparse_Matrix M     - the sparse matrix
 *      int row, col        - coordinates of element in matrix
 *      double AValue        - increment
 * Local variables:
 *      double *E        - pointer to the location in memory where the value
 *                        of this entry is stored
 * Algorithm:
 *      Call _SM_FindEntry to locate the memory address of the entry. It is
 *      called with _add_=TRUE (see above) so new entry will be allocated if
 *      necessary.
 *      Perform the operation using the pointer returned by _SM_FindEntry.
 ****************************************************************************/
{
    double *E;
    E=SM_FindEntry(M, row, col, TRUE);
    (*E)+=AValue;
    return;
}

double SM_GetEntry(Sparse_Matrix *M, int row, int col)
/****************************************************************************
 * Routine:     SM_GetEntry
 * ------------------------------------------------------------------------
 * Purpose:
 *      Returns the value of an entry in a sparse matrix
 * Parameters:
 *      Sparse_Matrix M     - the sparse matrix
 *      int row, col        - coordinates of element in matrix
 * Return value:
 *      double   the value of the element
 * Local variables:
 *      double *E        - pointer to the location in memory where the value
 *              of this entry is stored
 * Algorithm:
 *      Call _SM_FindEntry to locate the memory address of the entry. It is
 *      called with _add_=FALSE (see above) so if entry doesn't exist in memory
 *      (meaning the entry is 0) the returned value will be NULL.
 *      Return the value pointed to by the pointer returned by _SM_FindEntry or
 *      0.0 if pointer is NULL.
 ****************************************************************************/
{
    double *E=SM_FindEntry(M, row, col, FALSE);
    if(E)
        return *E;
    else
        return 0.0;
}



void MultiplySparseByVector(Sparse_Matrix *M, int dvect,
        const double *vect, int dres, double *res)
        /****************************************************************************
         * Routine:     MultiplySparseByVector
         * ------------------------------------------------------------------------
         * Purpose:
         *      Multiplies Sparse_Matrix by a vector
         * Parameters:
         *      Sparse_Matrix *M    - the matrix
         *      int dvect        - dimesion of the vector
         *      double *vect         - vector to be multiplied
         *      int dres         - dimesion of the resulting vector
         *      double *res          - the resulting vector
         * Return value: TRUE if successful, FALSE if dimestions disagree
         * Local variables:
         * Algorithm:
         *      Check correctness of input and call one of the above three routines
         *      to do the job
         ****************************************************************************/
{
    if((dvect!=M->cols)||(dres!=M->rows))
        error_msg("Cannot multiply matrix by vector.");
    if(M->data_type[0]==0)
    {
        if(M->Symmetric)
            _SM_MultSymmV_Sparse_Row((Sparse_Row*)M->SPARSE_DATA,
                    dres, vect, &res);
        else
            _SM_MultV_Sparse_Row((Sparse_Row*)M->SPARSE_DATA,
                    dres, vect, &res);
    }
    else
        _SM_MultV_Sparse_CSC(M->SPARSE_DATA, vect, &res);
    return;
}

void MultiplySparseTransByVector(Sparse_Matrix *M, int dvect,
        const double *vect, int dres, double *res)
        /****************************************************************************
         * Routine:     MultiplySparseTransByVector
         * ------------------------------------------------------------------------
         * Purpose:
         *      Multiplies the transpose of a Sparse_Matrix by a vector
         * Parameters:
         *      Sparse_Matrix *M    - the matrix
         *      int dvect        - dimesion of the vector
         *      double *vect         - vector to be multiplied
         *      int dres         - dimesion of the resulting vector
         *      double *res          - the resulting vector
         * Return value: TRUE if successful, FALSE if dimestions disagree
         * Local variables:
         * Algorithm:
         *      Check correctness of input and call one of the above three routines
         *      to do the job
         ****************************************************************************/
{
    if((dvect!=M->rows)||(dres!=M->cols))
        error_msg("Cannot multiply matrix by vector.");
    if(M->data_type[0]==0)
    {
        if(M->Symmetric)
            _SM_MultSymmV_Sparse_Row((Sparse_Row*)M->SPARSE_DATA,
                    dres, vect, &res);
        else
            _SM_MultTransV_Sparse_Row((Sparse_Row*)M->SPARSE_DATA,
                    dvect, dres, vect, &res);
    }
    else
        _SM_MultTransV_Sparse_CSC(M->SPARSE_DATA, vect, &res);
    return;
}

/*
 * int ForEachNonZeroEntryDo(Sparse_Matrix *M, SM_ENTRY_ACTION action,
 * void *CustomData)
 * {
 * int i, j;
 * Sparse_Row *R;
 * R = (Sparse_Row*)M->SPARSE_DATA;
 * if(M->Symmetric)
 * for(i=0;i<M->rows;i++, R++)
 * {
 * if(!action(i, R->Columns[0], R->Values[0], CustomData))
 * return FALSE;
 * for(j=1;j<R->EntriesCount;j++)
 * if((!action(i, R->Columns[j], R->Values[j], CustomData))
 * ||(!action(R->Columns[j], i, R->Values[j], CustomData)))
 * return FALSE;
 * }
 * else
 * for(i=0;i<M->rows;i++, R++)
 * for(j=0;j<R->EntriesCount;j++)
 * if(!action(i, R->Columns[j], R->Values[j], CustomData))
 * return FALSE;
 *
 * return TRUE;
 * }
 *
 *
 * void AddSparseMatrices(int n, double *C, Sparse_Matrix *M[], Sparse_Matrix
 *Res)
 * {
 * int i, r, c, s;
 * Sparse_Row *R;
 *
 * r=(*Res).rows;
 * c=(*Res).cols;
 * s=(*Res).Symmetric;
 * for(i=0;i<n;i++)
 * {
 * if((r!=M[i]->rows)||(c!=M[i]->cols))
 * error_msg("cannot add matrices.");
 * if(s&&(!M[i]->Symmetric))
 * error_msg("cannot add matrices.");
 * }
 * for(i=0;i<n;i++)
 * {
 * R = (Sparse_Row*)M[i]->SPARSE_DATA;
 * if((*Res).Symmetric==M[i]->Symmetric)
 * for(c=0;c<r;c++, R++)
 * for(s=0;s<R->EntriesCount;s++)
 * SM_AddToEntry(Res, c, R->Columns[s], C[i]*R->Values[s]);
 * else  // Res is not symmetric but M[i] is... //
 * for(c=0;c<r;c++, R++)
 * {
 * SM_AddToEntry(Res, c, c, C[i]*R->Values[0]);
 * for(s=1;s<R->EntriesCount;s++)
 * {
 * SM_AddToEntry(Res, c, R->Columns[s], C[i]*R->Values[s]);
 * SM_AddToEntry(Res, R->Columns[s], c, C[i]*R->Values[s]);
 * }
 * }
 * }
 * return;
 * }
 */
/*----------------------------------------------------------------------------*/

void PrintVector(double *V, int dim, char *filename)
{
    FILE *F;
    int i;
    if(V==NULL)
        return;
    if( (F=fopen(filename, "wt")) == NULL )
        return;
    for(i=0;i<dim;i++)
    {
        fprintf(F, "%-18.12g\n", (double) V[i]);
    }
    fclose(F);
    return;
}

void PrintSMatrix(Sparse_Matrix *M, char *filename)
{
    FILE *F;
    
    if( (M == NULL) )
        return;
    if ( (F = fopen(filename, "wt")) == NULL )
        return;
    fprintf(F, "%d,%d,%d\n",
            (int)M->rows, (int)M->cols, M->Symmetric);
    /* fprintf(F, "%s,%s,%s\n", "row", "col", "val"); */
    if(M->data_type[0]==0)
        PrintSRow(F, M->SPARSE_DATA, M->rows);
    else
        PrintSparse_CSC(F, M->SPARSE_DATA);
    fprintf(F, "\n");
    fclose(F);
}


/*----------------------------------------------------------------------------*/
/*  Linear systems */
/*----------------------------------------------------------------------------*/

/* Forward declarations */
int LUDecomposition(MATRIX M, int dim, MATRIX *LU);
int LUBackwardSubstitution(MATRIX LU, int dim, const double *rhs, VECTOR *s);
int LU_Solve(MATRIX M, int dim, const double *rhs, VECTOR *sol);
int Gauss_Solve(MATRIX _M, int dim, const double *_rhs, VECTOR *sol);
int GSI_Solve(MATRIX matr, int dim, const double *rhs, VECTOR *sol, double *TOL);

int GSI_Sparse(Sparse_Matrix *M, int dim, const double *rhs, VECTOR *sol, double
        *TOL);
void __sparse_solver_symm_mult(int dim, VECTOR A, VECTOR *B);
void __sparse_solver_mult(int dim, VECTOR A, VECTOR *B, int trans);
int CGM_Solve(void (*matrXvect)(int dim, VECTOR A, VECTOR *B),
        int dim, VECTOR rhs, VECTOR *sol,
        double TOL, double *RES, double *EPS);
int BCGM_Solve(void (*matrXvect)(int n, VECTOR A, VECTOR *B, int trans),
        int dim, VECTOR rhs, VECTOR *sol,
        double TOL, double *RES, double *EPS);


/*----------------------------------------------------------------------------*/
/*               LU-Decomposition (naive, i.e. no pivoting or scaling) */
/*----------------------------------------------------------------------------*/


int LUDecomposition(MATRIX M, int dim, MATRIX *LU)
/* dim - dimensions of the matrix */
/*  M - the matrix */
/*  LU - the LU-decomposed matrix */
{
    long i, j, k; /* counters */
    double tmp;  /* temporary variable */
    /*  first row of U is the same as the first row in M: */
    for(j=0; j<(long)dim; j++)
        (*LU)[0][j] = M[0][j];
    /*  rows loop starting from the second row */
    for(i=1; i<(long)dim; i++)
    {
        /*  fill L to the left of main diagonal */
        for(j=0; j<=i-1; j++)
        {
            for(tmp=.0, k=0; k<=j-1; k++)
            {
                tmp = tmp + (*LU)[i][k] * (*LU)[k][j];
            }
            (*LU)[i][j] = (M[i][j]-tmp)/(*LU)[j][j];
        }
        /*  Put U[i] on the main diagonal. */
        /* No need to store 1 on the main diagonal of L */
        for(tmp=.0, k=0; k<=i-1; k++)
            tmp = tmp + (*LU)[i][k] * (*LU)[k][i];
        /* if diagonal element is 0.0 then either pivoting is required, */
        /* but not implemented here, or the system is singular... */
        tmp = (*LU)[i][i] = M[i][i] - tmp;
        if (fabs(tmp)<SMALLF)
            return FALSE; /* indicates an error */
        /* fill U to the right of main diagonal */
        for(j=i+1; j<(long)dim; j++)
        {
            for(tmp=.0, k=0; k<=i-1; k++)
                tmp = tmp + (*LU)[i][k] * (*LU)[k][j];
            (*LU)[i][j] = M[i][j] - tmp;
        }
    }
    /*  if we are here then everything is fine */
    return TRUE; /* indicates success */
}


int LUBackwardSubstitution(MATRIX LU, int dim, const double *rhs, VECTOR *s)
/*  dim - dimension of linear system */
/*  LU - LU-decomposed nxn-dimensional matrix */
/*  rhs - n-dimensional righ hand side vector */
/*  s - n-dimensional solution vector */
{
    long i, j; /* counters */
    /* first solve L*y = rhs */
    /* s is used also as temporary variable "y" */
    for(i=0;i<(long)dim;i++)
    {
        (*s)[i] = rhs[i];
        for(j=0;j<i;j++)
            (*s)[i]-=LU[i][j]*(*s)[j];
    }
    /* then U*x = y */
    for(i=dim-1;i>=0;i--)
    {
        for(j=i+1;j<(long)dim;j++)
            (*s)[i]-=LU[i][j]*(*s)[j];
        (*s)[i] /= LU[i][i];
    }
    /* no problem can occur here provided lu_decomp was successful */
    return TRUE;
}

int LU_Solve(MATRIX M, int dim, const double *rhs, VECTOR *sol)
/*  dim - dimension of matrix */
/*  M - matrix of system */
/*  rhs - right hand side of system */
/*  sol - solution */
{
    MATRIX LU;          /* temporary storage for LU-decomposiotion */
    int result=FALSE;   /* return value of function -- FALSE by default */
    /* allocate memory for decomposed matrix - LU */
    LU = allocM(dim, dim);
    /* if LU decomposition returns FALSE then system is
     * singular or needs pivoting */
    if(!LUDecomposition(M, dim, &LU))
        goto clear_and_go;
    /* solve the system finally:
     * always returns 0 indicating success */
    result=LUBackwardSubstitution(LU, dim, rhs, sol);
    
    clear_and_go:
        /* release memory used for LU */
        freeM(&LU);
        return result;
}

/*----------------------------------------------------------------------------*/
/*               Gauss Direct Method */
/*----------------------------------------------------------------------------*/

int Gauss_Solve(MATRIX _M, int dim, const double *_rhs, VECTOR *sol)
{
    MATRIX M = allocM(dim, dim);
    double tmp, *v, *rhs = allocV(dim);
    long i, j, k;
    double *M0;
    
    M0 = M[0];
    for(i=0;i<(long)dim;i++)
        memcpy(M[i], _M[i], dim*sizeof(double));
    memcpy(rhs, _rhs, dim*sizeof(double));
    
    for(j=i=0;i<(long)dim-1;j=0, i++)
    {
        if(fabs(M[i][i]) < SMALLF) /* kind of pivoting... */
            for(j=i+1;j<(long)dim;j++)
                if(fabs(M[j][i]) > SMALLF)
                {
            v=M[i];
            M[i]=M[j];
            M[j]=v;
            tmp=rhs[i];
            rhs[i]=rhs[j];
            rhs[j]=tmp;
            break;
                }
        /* singular matrix -- pivoting didn't help */
        if(j==(long)dim)
            return FALSE; /* indicates error */
        
        for(j=i+1;j<(long)dim;j++)
        { /* Gaussian elimination: */
            tmp = M[j][i] / M[i][i];
            if(fabs(tmp) > SMALLF)
            {
                rhs[j] -= rhs[i] * tmp;
                for(k=i+1;k<(long)dim;k++)
                    M[j][k] -= M[i][k] * tmp;
            }
        }
    } /* Backward substitution */
    if(fabs(M[dim-1][dim-1]) < SMALLF)
        return FALSE;
    for(i=dim-1;i>=0;i--)
    {
        for(tmp=rhs[i], j=i+1; j<(long)dim; j++)
            tmp -= M[i][j]*(*sol)[j];
        (*sol)[i] = tmp / M[i][i];
    }
    /* claen before we leave: */
    M[0] = M0;
    freeM(&M);
    freeV(&rhs);
    return TRUE;
}

/*----------------------------------------------------------------------------*/
/*               Gauss-Seidel Iterative Method */
/*----------------------------------------------------------------------------*/

int GSI_Solve(MATRIX matr, int dim, const double *rhs, VECTOR *sol, double *TOL)
{
    const int MAX_ITER = 100;
    int iter, i, j;
    double inf_norm, sol_i;
    /*   double norm_2; */
    
    for(i=0;i<dim;i++)
        if(matr[i][i]==0.0)
            return FALSE;
    for(iter=0; iter<MAX_ITER; iter++)
    {
        inf_norm = 0.0;
        /*       norm_2 = 0.0; */
        for(i=0;i<dim;i++)
        {
            sol_i = rhs[i];
            for(j=0;j<dim;j++)
                sol_i -= matr[i][j] * (*sol)[j];
            sol_i = sol_i / matr[i][i];
            inf_norm = fmax(inf_norm, fabs(sol_i));
            /*           norm_2 += sol_i*sol_i; */
            (*sol)[i] += sol_i;
        }
        /*       norm_2 = sqrt(norm_2); */
        if(inf_norm<(*TOL))
            return iter;
    }
    (*TOL)=inf_norm;
    return FALSE;
}

int GSI_Sparse(Sparse_Matrix *M, int dim, const double *rhs, VECTOR *sol, double
        *TOL)
{
    const int MAX_ITER = 10000;
    int iter, i, j;
    double inf_norm, sol_i;
    /*   double norm_2; */
    int *Clmns;
    Sparse_Row *R, *Rows;
    R = Rows = (Sparse_Row*)M->SPARSE_DATA;
    
    for(i=0;i<dim;i++, R++)
        if(fabs(R->Values[0]) < SMALLF)
            error_msg("GSI_Sparse : zero on main diagonal.");
    for(iter=0; iter<MAX_ITER; iter++)
    {
        inf_norm = 0.0;
        /*       norm_2 = 0.0; */
        for(R=Rows, i=0;i<dim;i++, R++)
        {
            sol_i = rhs[i];
            Clmns = R->Columns;
            for(j=0;j<R->EntriesCount;j++)
                sol_i -= R->Values[j] * (*sol)[Clmns[j]];
            sol_i = sol_i / R->Values[0]; /*diagonal element is !=0 */
            inf_norm = fmax(inf_norm, fabs(sol_i));
            /*           norm_2 += sol_i*sol_i; */
            (*sol)[i] += sol_i;
        }
        /*       norm_2 = sqrt(norm_2); */
        if(inf_norm<(*TOL))
            return iter;
    }
    (*TOL)=inf_norm;
    return FALSE;
}
/*----------------------------------------------------------------------------*/
/*  Conjugate Gradient Method */
/*----------------------------------------------------------------------------*/
int CGM_Solve(void (*matrXvect)(int n, VECTOR A, VECTOR *B),
        int dim, VECTOR rhs, VECTOR *sol,
        double TOL, double *RES, double *EPS)
{
    static VECTOR r, p, tmp;
    static int vect_dim;
    double norm_p, alpha, beta, eps, n_sol, n_rhs;
    int i, iters;
    
    if(vect_dim < dim)
    {
        if(vect_dim)
        { freeV(&r); freeV(&p); freeV(&tmp); }
        r = allocV(dim);
        p = allocV(dim);
        tmp = allocV(dim);
        vect_dim = dim;
    }
    
    (*matrXvect)(dim, *sol, &tmp);
    for(n_rhs=(*RES)=i=0;i<dim;i++)
    {
        p[i]=r[i]=rhs[i]-tmp[i];
        n_rhs += rhs[i]*rhs[i];
        (*RES)+=p[i]*p[i];
    }
    n_rhs = sqrt(n_rhs);
    if(n_rhs <= TOL)
        n_rhs = 1.0;
    (*RES) = sqrt(*RES)/n_rhs;
    if((*RES)<=TOL)
    {
        (*EPS)=0.0;
        iters=0;
    }
    else
    {
        for(iters=1;;iters++)
        {
            (*matrXvect)(dim, p, &tmp);
            norm_p=n_sol=alpha=beta=(*RES)=(*EPS)=0.0;
            for(i=0;i<dim;i++)
            {
                norm_p+=p[i]*tmp[i];
                alpha+=r[i]*p[i];
            }
            alpha/=norm_p;
            for(i=0;i<dim;i++)
            {
                eps=alpha*p[i];
                (*sol)[i]+=eps;
                (*EPS)+=eps*eps;
                n_sol += (*sol)[i]*(*sol)[i];
                r[i]-=alpha*tmp[i];
                beta+=r[i]*tmp[i];
                (*RES)+=r[i]*r[i];
            }
            (*EPS)=sqrt(*EPS);
            if(n_sol>TOL)
                (*EPS)/=n_sol;
            (*RES)=sqrt(*RES)/n_rhs;
            for(beta/=norm_p, i=0;i<dim;i++)
                p[i]=r[i]-beta*p[i];
            if(((*EPS)<=TOL)||((*RES)<=TOL))
                break;
        }
    }
    return iters;
}


/*----------------------------------------------------------------------------*/
/*  Bi-Conjugate Gradient Method */
/*----------------------------------------------------------------------------*/
int BCGM_Solve(void (*matrXvect)(int n, VECTOR A, VECTOR *B, int trans),
        int dim, VECTOR rhs, VECTOR *sol,
        double TOL, double *RES, double *EPS)
{
    /* tmp will contain A.p */
    static VECTOR r, r_bar, p, p_bar, tmp;
    static int vect_dim;
    double p_dot, r_dot, alpha, beta, eps, norm_sol, n_rhs;
    int i, iters;
    
    if( vect_dim < dim )
    {
        freeV(&p);
        freeV(&p_bar);
        freeV(&r);
        freeV(&r_bar);
        freeV(&tmp);
        p=allocV(dim);
        p_bar=allocV(dim);
        r=allocV(dim);
        r_bar=allocV(dim);
        tmp=allocV(dim);
        vect_dim=dim;
    }
    
    (*matrXvect)(dim, *sol, &tmp, FALSE);
    for(n_rhs=r_dot=i=0;i<dim;i++)
    {
        p_bar[i]=r_bar[i]=p[i]=r[i]=rhs[i]-tmp[i];
        r_dot += r_bar[i] * r[i]; /* same as r[i] * r[i] ! */
        n_rhs += rhs[i] * rhs[i];
    }
    n_rhs = sqrt(n_rhs);
    if(n_rhs < TOL)
        n_rhs = 1.0;
    (*RES) = sqrt(r_dot) / n_rhs;
    if((*RES)<TOL)
    {
        (*EPS)=0.0;
        iters=0;
    }
    else
    {
        for(iters=1;iters<=dim*10;iters++)
        {
            (*matrXvect)(dim, p, &tmp, FALSE);
            (*RES)=(*EPS)=norm_sol=p_dot=0.0;
            for(i=0;i<dim;i++)
                p_dot+=p_bar[i]*tmp[i];
            alpha=r_dot/p_dot;
            for(i=0;i<dim;i++)
            {
                eps=alpha*p[i];
                (*EPS) += eps*eps;
                (*sol)[i] += eps;
                norm_sol += (*sol)[i] * (*sol)[i];
                r[i] -= alpha * tmp[i];
                (*RES) += r[i] * r[i];
            }
            norm_sol = sqrt(norm_sol);
            (*EPS)=sqrt(*EPS);
            if(norm_sol>TOL)
                (*EPS)/=norm_sol;
            (*RES)=sqrt(*RES)/n_rhs;
            if(((*EPS)<=TOL)||((*RES)<=TOL))
                break;
            (*matrXvect)(dim, p_bar, &tmp, TRUE);
            for(i=0;i<dim;i++)
                r_bar[i] -= alpha * tmp[i];
            beta = 1.0 / r_dot;
            for(r_dot=i=0;i<dim;i++)
                r_dot += r_bar[i] * r[i];
            beta *= r_dot;
            for(i=0;i<dim;i++)
            {
                p[i] = r[i] + beta * p[i];
                p_bar[i] = r_bar[i] + beta * p_bar[i];
            }
        }
    }
    return iters;
}


/*----------------------------------------------------------------------------*/
/*               Linear System Solvers */
/*----------------------------------------------------------------------------*/

int SolveLinearSystem(const MATRIX M, int dim, const double *rhs, VECTOR sol)
{
    double tol = SOLVER_TOLERANCE;
    int res;
    res = Gauss_Solve(M, dim, rhs, &sol);
    if(res)
        return res;
    res = LU_Solve(M, dim, rhs, &sol);
    if(res)
        return res;
    memset(sol, 0, dim*sizeof(double));
    if((res=GSI_Solve(M, dim , rhs, &sol, &tol))!=0)
        return res;
    return res;
}


/*Sparse_Row *__sparse_solver_matrix_row_data=NULL;*/
/* this one is used for the CGM */
/*void __sparse_solver_symm_mult(int dim, VECTOR A, VECTOR *B)
 * {
 * if(__sparse_solver_matrix_row_data)
 * _SM_MultSymmV(__sparse_solver_matrix_row_data,dim,A,B);
 * return;
 * }*/
/* this one is for non-symmetric matrices */
/*void __sparse_solver_mult(int dim, VECTOR A, VECTOR *B, int trans)
 * {
 * if(__sparse_solver_matrix_row_data)
 * {
 * if( trans )
 * _SM_MultTransV(__sparse_solver_matrix_row_data,dim,dim,A,B);
 * else
 * _SM_MultV(__sparse_solver_matrix_row_data,dim,A,B);
 * }
 * return;
 * }*/

#include <time.h>
/*
 * int SolveSparseSystem(Sparse_Matrix *M, int dim, VECTOR rhs, VECTOR *sol,
 * double *RES, double *ERR)
 * {
 * int iters;
 * double res, eps;
 * clock_t tm;
 *
 * PrintLog("CGM: ");
 * if(M->rows!=M->cols)
 * error_msg("Matrix must be square.");
 * if(dim!=M->rows)
 * error_msg("Dimensions of vectors and matrix disagree.");
 * PrintLog(" dim=%-7u  ", M->rows);
 * PrintLogdbg(("    Target accuracy    %g\n", SOLVER_TOLERANCE));
 * tm = clock();
 * if(!M->Symmetric)
 * {
 * //
 * //        eps = SOLVER_TOLERANCE; res = 0.0;
 * //        iters=GSI_Sparse(M, dim, rhs, sol, &eps);
 * //
 * __sparse_solver_matrix_row_data = (Sparse_Row*)M->SPARSE_DATA;
 * iters=BCGM_Solve(__sparse_solver_mult, dim, rhs, sol,
 * SOLVER_TOLERANCE, &res, &eps);
 * __sparse_solver_matrix_row_data = NULL;
 * }
 * else
 * {
 * __sparse_solver_matrix_row_data = (Sparse_Row*)M->SPARSE_DATA;
 * iters=CGM_Solve(__sparse_solver_symm_mult, dim, rhs, sol,
 * SOLVER_TOLERANCE, &res, &eps);
 * __sparse_solver_matrix_row_data = NULL;
 * }
 * tm = clock() - tm;
 * if(RES)
 * (*RES)=res;
 * if(ERR)
 * (*ERR)=eps;
 * PrintLog("iters=%-6u  eps=%-8.2e  res=%-8.2e  CPU time=%-4.2fsec.\n",
 * iters, (double) eps, (double) res, tm / (double) CLOCKS_PER_SEC);
 * return iters;
 * }
 */

/*--------------------------------------------------------------------------*/
/*  Conjugate Gradient Method */
/*--------------------------------------------------------------------------*/
static void ps_null(double *X, const double *F, void *ps_data)
{
    (void) X; (void)F; (void)ps_data;
    return;
}

/*************************   cgm   ******************************************/
void solve_cgm(matr_prod_t matr_prod, prec_solv_t prec_solv,
        int dim, double *X, const double *F, void *mp_data, void *ps_data,
        double TOL, double *RES, double *EPS, int *iters,
        double *CPU)
{
    double *R = NULL, *P, *T, *Z;
    double nP, alpha, beta, delta, nX, nF;
    double res, eps, cpu;
    int i, it;
    clock_t tm0, tm1;
    
    cpu = 0.0;
    tm0 = clock();
    
    if(dim==0)
        error_msg("Matrix of size 0x0 passed to 'solve_cgm'\n");
    
    /* compute the norm of the right-hand side */
    for(nF=i=0;i<dim;i++)
        nF += F[i] * F[i];
    nF = sqrt(nF);
    
    if(nF <= SMALLF) /* trivial solution */
    {
        CleanVector(X, dim);
        eps = res = it = 0;
        goto leave;
    }
    
    TOL = TOL<SMALLF ? SOLVER_TOLERANCE : TOL;
    nF = fmax(TOL, nF);
    /* allocate memory.  after this do not use 'return' as we release
     * the memory at the end. */
    if(prec_solv)
        Z = (R = allocV(4*dim)) + dim;
    else
    {
        Z = R = allocV(3*dim);
        prec_solv = ps_null;
    }
    P = Z + dim;
    T = P + dim;
    
    matr_prod(X, T, FALSE, mp_data);
    for(res=i=0;i<dim;i++)
    {
        P[i]=R[i]=F[i]-T[i];
        res+=R[i]*R[i];
    }
    res = sqrt(res)/nF;
    if(res <= SMALLF)
    {             /* initial guess is the solution, yay */
        eps = it = 0;
        goto leave;
    }
    for(it=1;;it++)
    {
        matr_prod(P, T, FALSE, mp_data);
        nP=nX=alpha=beta=res=eps=0.0;
        for(i=0;i<dim;i++)
        {
            nP+=P[i]*T[i];
            alpha+=R[i]*P[i];
        }
        alpha /= nP;
        for(i=0;i<dim;i++)
        {
            delta=alpha*P[i];
            X[i]+=delta;
            eps+=delta*delta;
            nX += sqr(X[i]);
            R[i]-=alpha*T[i];
            res+=R[i]*R[i];
        }
        nX = fmax(TOL, nX);
        eps=sqrt(eps) / nX;
        res=sqrt(res) / nF;
        if((res <= SMALLF) || (eps <= SMALLF) || (res <= TOL && eps <= TOL))
            /* if((res <= TOL && eps <= TOL)) */
            break;
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        tm0 = tm1;
        prec_solv(Z, R, ps_data);
        for(i=0;i<dim;i++)
        {
            beta+=Z[i]*T[i];
        }
        beta /= nP;
        for(i=0;i<dim;i++)
            P[i]=Z[i]-beta*P[i];
    }
    leave:
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        if(CPU) *CPU = cpu;
        if(EPS) *EPS = eps;
        if(RES) *RES = res;
        if(iters) *iters = it;
        if(R) freeV(&R);
}

void __cgm_by_shewchuk_94(matr_prod_t matr_prod, prec_solv_t prec_solv,
        int dim, double *X, const double *F, void *mp_data, void *ps_data,
        double TOL, double *RES, double *EPS, int *iters,
        double *CPU)
{
    double *R, *P, *T, *Z;
    double D0, DN, DO, alpha, beta;
    int i, RRI;
    double nF, nX;
    double cpu;
    clock_t tm0, tm1;
    
    cpu = 0.0;
    tm0 = clock();
    
    RRI = max(50, floor(sqrt(dim)));
    TOL = TOL<SMALLF ? SOLVER_TOLERANCE : TOL;
    
    /* allocate memory first */
    if(prec_solv==NULL)
        Z = R = (double *)calloc(3*dim, sizeof(double));
    else
        Z = (R = (double *)calloc(4*dim, sizeof(double))) + dim;
    if(R==NULL)
        error_msg("solve_cgm: memory allocation error.\n");
    P = Z + dim;
    T = P + dim;
    
    /* norm of F */
    nF = 0.0;
    for(i=0;i<dim;i++)
        nF += F[i]*F[i];
    nF = sqrt(nF);
    if(nF < TOL)
        nF = 1.0;
    
    *iters = 0;
    /* R = F - A*X */
    matr_prod(X, T, FALSE, mp_data);
    for(i=0;i<dim;i++)
        R[i] = F[i] - T[i];
    /* P = B^(-1) * R */
    if(prec_solv)
        prec_solv(P, R, ps_data);
    else
        memcpy(P, R, dim * sizeof(double));
    /* DN = R * P */
    DN = 0.0;
    for(i=0;i<dim;i++)
        DN += R[i] * P[i];
    /* */
    D0 = DN;
    while(DN > D0 * TOL * TOL)
    {
        /* T = A * P */
        matr_prod(P, T, FALSE, mp_data);
        /* alpha = DN / (P * T) */
        alpha = 0.0;
        for(i=0;i<dim;i++)
            alpha += P[i] * T[i];
        alpha = DN / alpha;
        /* X = X + alpha * P */
        for(nX=0.0, i=0;i<dim;i++)
            nX += sqr(X[i] += alpha * P[i]);
        nX = sqrt(nX);
        if(nX <= TOL)
            nX = 1.0;
        /* update R */
        (*iters)++;
        if( ((*iters) % RRI) == 0 )
        {
            /* R = F - A*X */
            matr_prod(X, R, FALSE, mp_data);
            for(i=0;i<dim;i++)
                R[i] = F[i] - R[i];
        }
        else
        {
            /* R = R - alpha * T */
            for(i=0;i<dim;i++)
                R[i] -= alpha * T[i];
        }
        /* Z = B^(-1) * R */
        if(prec_solv)
            prec_solv(Z, R, ps_data);
        /* DO = DN */
        DO = DN;
        /*DN = R * Z */
        DN = 0.0;
        for(i=0;i<dim;i++)
            DN += R[i] * Z[i];
        /* beta = DN / DO */
        beta = DN / DO;
        /* compute norm of residual in RES */
        *RES = 0.0;
        for(i=0;i<dim;i++)
            *RES += sqr(R[i]);
        *RES = sqrt(*RES) / nF;
        /* compute norm of correction in EPS */
        *EPS = 0.0;
        for(i=0;i<dim;i++)
            *EPS += sqr(alpha * P[i]);
        *EPS = sqrt(*EPS) / nX;
        if((*EPS < TOL) || (*RES < TOL))
            break;
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        tm0 = tm1;
        /* P = Z + beta * P */
        for(i=0;i<dim;i++)
            P[i] = Z[i] + beta * P[i];
    }
    free(R);
    tm1 = clock();
    cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
    if(CPU) *CPU = cpu;
}

/****************************  bcgm  ****************************************/

void solve_bcgm(matr_prod_t matr_prod, prec_solv_t prec_solv,
        int dim, double *X, const double *F, void *mp_data, void *ps_data,
        double TOL, double *RES, double *EPS, int *iters,
        double *CPU)
{
    /* const double BREAKDOWN = 1e-9; */
    double *R = NULL, *Rb, *P, *Pb, *T;
    double nP, nR, nRold, alpha, beta, delta, nX, nF;
    double res, eps, cpu;
    int i, it;
    clock_t tm0, tm1;
    
    (void) ps_data;
    if(prec_solv)
        warning_msg("preconditioning not available in bcgm.  Using "
                "standard method.\n");
    
    if(dim==0)
        error_msg("Matrix of size 0x0 passed to 'solve_bcgm'\n");
    
    cpu = 0.0;
    tm0 = clock();
    
    eps = res = it = 0;
    /* compute the norm of the right-hand side */
    for(nF=i=0;i<dim;i++)
        nF += F[i] * F[i];
    nF = sqrt(nF);
    
    if(nF <= SMALLF) /* trivial solution */
    {
        CleanVector(X, dim);
        goto leave;
    }
    
    TOL = TOL<SMALLF ? SOLVER_TOLERANCE : TOL;
    nF = fmax(TOL, nF);
    
    /* allocate memory.  after this do not use 'return' as we release
     * the memory at the end. */
    R = allocV(5*dim);
    Rb = R + dim;
    P = Rb + dim;
    Pb = P + dim;
    T = Pb + dim;
    
    matr_prod(X, T, FALSE, mp_data);
    for(nR=i=0;i<dim;i++)
    {
        Pb[i]=P[i]=Rb[i]=R[i]=F[i]-T[i];
        nR+=Rb[i]*R[i];
    }
    /* nR = fmax(TOL, nR); */
    res = sqrt(nR) / nF;
    if(res <= SMALLF)
        goto leave;     /* initial guess is the solution, yay */
    
    for(it=1;it<dim*1.5;it++)
    {
        matr_prod(P, T, FALSE, mp_data);
        nRold=nR;
        nR=nP=nX=res=eps=0.0;
        for(i=0;i<dim;i++)
            nP+=Pb[i]*T[i];
        /*
         * if (fabs(nP) < BREAKDOWN * fabs(nRold))
         * error_msg("bcgm: breakdown at iteration %d. "
         * "nP=%-.4e is too small.\n", it, nP);
         */
        alpha = nRold / nP;
        for(i=0;i<dim;i++)
        {
            delta=alpha*P[i];
            X[i]+=delta;
            eps+=delta*delta;
            nX += sqr(X[i]);
            R[i]-=alpha*T[i];
            res+=R[i]*R[i];
        }
        nX = fmax(TOL, nX);
        eps=sqrt(eps) / nX;
        res=sqrt(res) / nF;
        if((res <= SMALLF) || (eps <= SMALLF) || (res <= TOL && eps <= TOL))
            /* if((res <= TOL && eps <= TOL)) */
            break;
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        tm0 = tm1;
        matr_prod(Pb, T, TRUE, mp_data);
        for(i=0;i<dim;i++)
        {
            Rb[i]-=alpha*T[i];
            nR+=Rb[i]*R[i];
        }
        /*
         * if (fabs(nRold) < BREAKDOWN * fabs(nR))
         * error_msg("bcgm: breakdown at iteration %d. "
         * "nRold=%-.4e is too small.\n", it, nRold);
         */
        beta = nR / nRold;
        for(i=0;i<dim;i++)
        {
            P[i]=R[i]+beta*P[i];
            Pb[i]=Rb[i]+beta*Pb[i];
        }
    }
    leave:
        if(it >= 1.5*dim)
            warning_msg("maximum number of iterations exceeded.\n");
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        if(CPU) *CPU = cpu;
        if(EPS) *EPS = eps;
        if(RES) *RES = res;
        if(iters) *iters = it;
        if(R) freeV(&R);
}


void solve_bcgstab(matr_prod_t matr_prod, prec_solv_t prec_solv,
        int dim, double *X, const double *F, void *mp_data, void *ps_data,
        double TOL, double *RES, double *EPS, int *iters,
        double *CPU)
        /*
         * This implementation of Bi-CGSTAB is taken from sec. 3.6.3 of
         *    C.T.Kelley, Iterative methods for linear and nonlinear equations. 1995
         */
{
    const double BREAKDOWN = SMALLF;
    double *R = NULL, *R0, *P, *V, *T, *S;
    double nR, alpha, beta, delta, nX, nF, omega, tau, ttau, nT, rho0, rho1;
    double res, eps, cpu;
    int i, it;
    clock_t tm0, tm1;
    
    (void) ps_data;
    if(prec_solv)
        warning_msg("preconditioning not available in bcgstab.  Using "
                "standard method.\n");
    
    if(dim==0)
        error_msg("Matrix of size 0x0 passed to '%s'\n", __func__);
    
    cpu = 0.0;
    tm0 = clock();
    
    eps = res = it = 0;
    /* compute the norm of the right-hand side */
    for(nF=i=0;i<dim;i++)
        nF += F[i] * F[i];
    nF = sqrt(nF);
    
    if(nF <= SMALLF * dim) /* trivial solution */
    {
        CleanVector(X, dim);
        goto leave;
    }
    
    TOL = TOL<SMALLF ? SOLVER_TOLERANCE : TOL;
    
    /* allocate memory.  after this do not use 'return' as we release
     * the memory at the end. */
    R = allocV(5*dim);
    P = R + dim;
    V = P + dim;
    R0 = V + dim;
    T = R0 + dim;
    S = R;
    
    matr_prod(X, T, FALSE, mp_data);
    nR = 0.0;
    for(i=0;i<dim;i++)
    {
        R0[i] = R[i] = F[i] - T[i];
        V[i] = P[i] = 0.0;
        nR += R0[i] * R[i];
    }
    rho0 = alpha = omega = 1.0;
    rho1 = nR;
    
    
    res = sqrt(nR) / fmax(TOL, nF);
    if(res < SMALLF)
        goto leave;     /* initial guess is the solution, yay */
    
    for(it=1;it<dim*1.5;it++)
    {
        if(fabs(omega) < BREAKDOWN * fabs(alpha))
            error_msg("bcgs: breakdown at iteration %d. "
                    "omega=%-.4e is too small.\n", it, omega);
        /* step 2b */
        beta = (rho1 / rho0) * (alpha / omega);
        /* step 2c */
        for(i=0;i<dim;i++) P[i] = R[i] + beta*(P[i] - omega*V[i]);
        /* step 2d */
        matr_prod(P, V, FALSE, mp_data);
        /* step 2e */
        tau = 0.0;
        for(i=0;i<dim;i++) tau += R0[i] * V[i];
        if(fabs(tau) < BREAKDOWN * fabs(rho1))
            error_msg("bcgs: breakdown at iteration %d. "
                    "tau=%-.4e is too small.\n", it, tau);
        alpha = rho1 / tau;
        /* step 2f */
        for(i=0;i<dim;i++) S[i] = R[i] - alpha*V[i];
        matr_prod(S, T, FALSE, mp_data);
        /* step 2g */
        nT = 0.0; tau = 0.0; ttau = 0.0;
        for(i=0;i<dim;i++)
        {
            nT += T[i] * T[i];
            tau += T[i] * S[i];
            ttau += R0[i] * T[i];
        }
        if( nT < BREAKDOWN * fabs(tau))
            error_msg("bcgs: breakdown at iteration %d. "
                    "nT=%-.4e is too small.\n", it, nT);
        omega = tau / nT;
        rho0 = rho1;
        rho1 = -omega * ttau;
        /* steps 2h and 2i*/
        nX = eps = nR = 0.0;
        for(i=0;i<dim;i++)
        {
            /* step 2h */
            delta = alpha*P[i] + omega*S[i];
            eps += delta*delta;
            X[i] = X[i] + delta;
            nX += X[i]*X[i];
            /* step 2i */
            R[i] = S[i] - omega * T[i];
            nR += R[i] * R[i];
        }
        /*
         * eps=sqrt(eps) / fmax(TOL, nX);
         * res=sqrt(nR) / fmax(TOL, nF);
         */
        /* the termination condition is different from the one in the book.
         * It is instead made similar to termination conditions in the other
         * iterative solvers here */
        eps=sqrt(eps);
        res=sqrt(nR);
        if((res < TOL && eps < TOL))
            break;
        
        /*
         * if(sqrt(nR) < TOL * nF)
         * break;
         */
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        tm0 = tm1;
    }
    leave:
        if(it >= 1.5*dim)
            warning_msg("maximum number of iterations exceeded.\n");
        tm1 = clock();
        cpu += (tm1-tm0) / (double) CLOCKS_PER_SEC;
        if(CPU) *CPU = cpu;
        if(EPS) *EPS = eps;
        if(RES) *RES = res;
        if(iters) *iters = it;
        if(R) freeV(&R);
}


/****************************************************************************/
void __mp_sparse(const double *X, double *F, int trans, void *mp_data)
{
    Sparse_Matrix *M = mp_data;
    if(trans)
        MultiplySparseTransByVector(M, M->cols, X, M->rows, F);
    else
        MultiplySparseByVector(M, M->cols, X, M->rows, F);
    
    /*
     * if(M->Symmetric)
     * _SM_MultSymmV((Sparse_Row*)M->SPARSE_DATA, M->rows, X, &F);
     * else
     * {
     * if(trans)
     * _SM_MultTransV((Sparse_Row*)M->SPARSE_DATA, M->rows, M->cols, X, &F);
     * else
     * _SM_MultV((Sparse_Row*)M->SPARSE_DATA, M->rows, X, &F);
     * }
     */
}

typedef struct
{
    int dim;
    double *D;
} __ps_diag_data_t;


void __ps_diag(double *X, const double *F, void *ps_data)
{
    __ps_diag_data_t *data = ps_data;
    int i;
    for(i=0;i<data->dim;i++)
        X[i] = F[i] / data->D[i];
}

void sparse_cgm(Sparse_Matrix *M, double *X, const double *F, double TOL)
{
    double R, E;
    int I;
    double CPU;
    
    if(M->rows != M->cols)
        error_msg("sparse_cgm: matrix must be square.\n");
    if(M->Symmetric)
    {
        PrintLog(" cgm:  dim=%-7u  ", M->rows);
        solve_cgm(__mp_sparse, NULL, M->rows, X, F, M, NULL, TOL,
                &R, &E, &I, &CPU);
    }
    else
    {
        /*
         * int i,j;
         * MATRIX full;
         * clock_t tm0;
         * Sparse_Row *Row;
         *
         * PrintLog("FULL:  dim=%-7u  ", M->rows);
         * tm0 = clock();
         * full = allocM(M->rows,M->cols);
         *
         * Row = (Sparse_Row*)M->SPARSE_DATA;
         * for(i=0;i<M->rows;i++,Row++)
         * for(j=0;j<Row->EntriesCount;j++)
         * full[i][Row->Columns[j]] = Row->Values[j];
         *
         * SolveLinearSystem(full,M->rows,F,X);
         * I = E = R = 0.0;
         *
         * freeM(&full);
         * CPU = (clock()-tm0) / (double) CLOCKS_PER_SEC;
         */
        
        
        /*
         * clock_t tm0;
         * PrintLog(" GSI:  dim=%-7u  ", M->rows);
         * E = TOL<SMALLF ? SOLVER_TOLERANCE : TOL;;
         * R = 0.0;
         * tm0 = clock();
         * I=GSI_Sparse(M, M->rows, F, &X, &E);
         * CPU = (clock()-tm0) / (double) CLOCKS_PER_SEC;
         */
        
        PrintLog("bcgm:  dim=%-7u  ", M->rows);
        solve_bcgm(__mp_sparse, NULL, M->rows, X, F, M, NULL, TOL,
                &R, &E, &I, &CPU);
        /*
         *
         * PrintLog("bcgs:  dim=%-7u  ", M->rows);
         * solve_bcgstab(__mp_sparse, NULL, M->rows, X, F, M, NULL, TOL,
         * &R, &E, &I, &CPU);
         */
        
    }
    PrintLog("iters=%-6u  eps=%-8.2e  res=%-8.2e  CPU time=%-4.2fsec.\n",
            I, (double) E, (double) R, CPU);
}


void sparse_pcgm_diag(Sparse_Matrix *M, double *P, double *X, const double *F,
        double TOL)
{
    __ps_diag_data_t data;
    double R, E;
    int I;
    double CPU;
    
    if(M->rows != M->cols || !M->Symmetric)
        error_msg("sparse_cgm: matrix must be square and symmetric.\n");
    PrintLog("pcgm:  dim=%-7u  ", M->rows);
    data.dim = M->rows;
    data.D = P;
    solve_cgm(__mp_sparse, __ps_diag, M->rows, X, F, M, &data, TOL,
            &R, &E, &I, &CPU);
    PrintLog("iters=%-6u  eps=%-8.2e  res=%-8.2e  CPU time=%-4.2fsec.\n",
            I, (double) E, (double) R, CPU);
}

void sparse_bcgm(Sparse_Matrix *M, double *X, const double *F, double TOL)
{
    double R, E;
    int I;
    double CPU;
    
    if(M->rows != M->cols)
        error_msg("sparse_bcgm: matrix must be square.\n");
    PrintLog("bcgm:  dim=%-7u  ", M->rows);
    solve_bcgm(__mp_sparse, NULL, M->rows, X, F, M, NULL, TOL,
            &R, &E, &I, &CPU);
    PrintLog("iters=%-6u  eps=%-8.2e  res=%-8.2e  CPU time=%-4.2fsec.\n",
            I, (double) E, (double) R, CPU);
}

static void SolveCSC(Sparse_CSC *data, double *X,
        const double *F, double TOL)
{
    #if !defined(USE_SPSOLVER)
        error_msg("%s: USE_SPSOLVER is not defined.\n", __func__);
    #else
    /*
     long neqns = data->neqns;
     long   nnz = data->colstr[neqns];
     
     long msglvl = 0;
     long nrhs = 1;
     long ldrhs = neqns;
     long outunt = 5;
     long ierr = 0;
     char pivot[] = "n";
     
     long i;
     
     */
    
    CopyVector(F, X, data->neqns);
    
    
    /*
     for(i=0;i<=neqns;i++) data->colstr[i]++;
     for(i=0;i<nnz;i++) data->rowind[i]++;
     */
    long one = 1, ierr = 0;
    assert(data->handle!=NULL);
    dgsssl_64_(&one, X, &data->neqns, data->handle, &ierr);
    if(ierr != 0) error_msg("%s: ierr = %d", "dgsssl_64_", ierr);
    
    
    /*
     assert(data->handle==NULL);
     data->handle = allocV(150);
     
     PrintLog("calling %s\n", "dgssin_64_");
     dgssin_64_(data->mtrtyp, pivot, &neqns, data->colstr, data->rowind,
     &outunt, &msglvl, data->handle, &ierr, 2, 1);
     if(ierr != 0) error_msg("%s: ierr = %d", "dgssin_64_", ierr);
     
     PrintLog("calling %s\n", "dgssor_64_");
     dgssor_64_(ordmth, data->handle, &ierr, 3);
     if(ierr != 0) error_msg("%s: ierr = %d", "dgssor_64_", ierr);
     
     PrintLog("calling %s\n", "dgssfa_64_");
     dgssfa_64_(&neqns, data->colstr, data->rowind, data->values,
     data->handle, &ierr);
     if(ierr != 0) error_msg("%s: ierr = %d", "dgssfa_64_", ierr);
     
     long one = 1;
     PrintLog("calling %s\n", "dgsssl_64_");
     dgsssl_64_(&one, X, &neqns, data->handle, &ierr);
     if(ierr != 0) error_msg("%s: ierr = %d", "dgsssl_64_", ierr);
     
     PrintLog("calling %s\n", "dgssda_64_");
     dgssda_64_(data->handle, &ierr);
     if(ierr != 0) error_msg("%s: ierr = %d", "dgssda_64_", ierr);
     
     free(data->handle);
     data->handle = NULL;
     */
    
    /*
     for(i=0;i<=neqns;i++) data->colstr[i]--;
     for(i=0;i<nnz;i++) data->rowind[i]--;
     */
    #endif
}

void SolveSparse(Sparse_Matrix *M, double *X, const double *F, double TOL)
{
    if(M->data_type[0]==0)
    {
        sparse_cgm(M, X, F, TOL);
        return;
    }
    else
    {
        #if defined(USE_SPSOLVER)
        SolveCSC(M->SPARSE_DATA, X, F, TOL);
        #else
        sparse_cgm(M, X, F, TOL);
        #endif
    }
}
